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ABSTRACT 

We analyze the global hydrodynamic flow in the ocean of an accreting, rapidly rotating, non-magnetic neutron 
star in a low-mass X-ray binary during a type I X-ray burst. We use both analytical arguments and numerical 
simulations of simplified models for ocean burning. Our analysis extends previous work by taking into account 
the rapid rotation of the star and the lift-up of the burning ocean during the burst. We find a new regime for the 
spreading of a nuclear burning front, where the flame is carried along a coherent shear flow across the front. If 
turbulent viscosity is weak, the speed of flame propagation is Vfl am e ~ (gh) 1 ' 2 //^ ~ 20 km s" 1 , where h is the 
scale height of the burning ocean, g is the local gravitational acceleration, t„ is the timescale for fast nuclear 
burning during the burst, and / is the Coriolis parameter, i.e., twice the local vertical component of the spin 
vector. If turbulent viscosity is dynamically important, the flame speed increases, and reaches the maximum value, 
v flame ~ (<?V/ f «) I//2 ~ 300 km s" 1 , when the eddy overturn frequency is comparable to the Coriolis parameter/. 

We show that, due to rotationally reduced gravity, the thermonuclear runaway which ignites the ocean is likely 
to begin on the equator. The equatorial belt is ignited at the beginning of the burst, and the flame then propagates 
from the equator to the poles. Inhomogeneous cooling (equator first, poles second) of the hot ashes drives strong 
zonal currents which may be unstable to the formation of Jupiter-type vortices; we conjecture that these vortices 
are responsible for coherent modulation of X-ray flux in the tails of some bursts. We consider the effect of strong 
zonal currents on the frequency of modulation of the X-ray flux and show that the large values of the frequency 
drifts observed in some bursts can be accounted for within our model combined with the model of homogeneous 
radial expansion. Additionally, if vortices or other inhomogeneities are trapped in the forward zonal flows around 
the propagating burning front, fast chirps with large frequency ranges (~ 25 - 500 Hz) may be detectable during 
the burst rise. Finally, we argue that an MHD dynamo within the burning front can generate a small-scale magnetic 
field, which may enforce vertically rigid flow in the front's wake and can explain the coherence of oscillations in 
the burst tail. 

Subject headings: accretion, accretion disks - instabilities - hydrodynamics - nuclear reactions, nucleosynthesis, 
abundances - stars: neutron - stars: rotation - X-rays: bursts 



1. INTRODUCTION 

Accreting neutron stars (NSs) in low-mass X-ray bina- 
ries (LMXBs) undergo type I X-ray bursts due to thermonu- 
clear runaways in pure helium or mixed hydrogen/helium lay- 
ers (Hansen and Van Horn 1975, Maraschi and Cavaliere 
1976; Woosley and Taam 1996; see Lewin, van Paradijs, and 
Taam 1995 and Bildsten 1998 for reviews). Spherically sym- 
metric models of such runaways successfully explain gen- 
eral features of type I X-ray bursts, such as burst fluences 
(~ 10 39 erg), timescales for accretion between bursts (~ few 
hours), and burst durations (~ 10 seconds). However, spheri- 
cally symmetric models cannot account for the lateral spreading 
of thermonuclear flames and its interplay with NS rotation, and 
recent observations have brought this issue into focus. 

Since the timescale for accretion between bursts is much 
longer than the burst duration, it is unlikely that identical con- 
ditions exist over the whole stellar surface for the burning insta- 
bility to start simultaneously (Shara 1981). Therefore, burning 
should start locally at some point, creating a brightness asym- 
metry, and spread over the entire surface of the star. As nuclear 
burning spreads around the star, one would then expect to see a 
rotational modulation of the X-ray flux, with the frequency of 



the modulation equal that of the neutron-star spin. Such highly 
coherent (Q ~ few thousand) "burst oscillations" have indeed 
been observed with RXTE in nine different bursters for many 
bursts with rise times of < lsec (see van der Klis 2000 for a 
review and Muno et. al. 2001 for the most recent tally). The 
inferred neutron-star spins v s are between 250 and 650 Hz. 

The discovery of burst oscillations, while confirming the 
basic expectations of millisecond spins of accreting NSs and 
asymmetry of nuclear burning, has brought about a host of 
new questions. Firstly, burst oscillations are seen only from 
some of the ~ 50 known Galactic LMXBs, and only some 
bursts from the same source show oscillations. For NSs with 
v s ~ 600 Hz, oscillations are seen only during strong bursts 
with photospheric radius expansion, while for v s ~ 300 Hz, os- 
cillations are equally as likely to be seen during weak or strong 
bursts (Muno et. al. 2001). Secondly, oscillations are most com- 
monly seen during tails of bursts, when, presumably, the entire 
accreted fuel layer has been burned and the obvious asymmetry 
is no longer present. Finally, the frequency of burst oscillations 
drifts upward by Av ~ several Hz during the burst. 

A simple model for the drift of the burst oscillation frequency 
has been proposed by Strohmayer et al. (1997), and recently 
considered quantitatively by Cumming and Bildsten (2000) and 
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Cumming et. al. (2001). Since the vertical sound crossing time 
through the burning layer (microseconds) is much smaller than 
the nuclear burning time (> fraction of a second) on which 
the burst evolves, the outer layers of the NS are always in hy- 
drostatic balance. When the burning layer is hot, it expands 
hydrostatically by Az ^tens of meters while conserving an- 
gular momentum, and, hence, lags behind the neutron star by 
Av/v s w 2Az/R. During the tail of the burst, a postulated 
temperature inhomogeneity gives rise to oscillations, and, as 
the layer cools down and contracts, we observe an upward fre- 
quency drift by Av of a few Hz, roughly consistent with obser- 
vations. However, recent works by Van Straaten et al. (2001), 
Galloway et al. (2001), Wijnands et al. (2001), and Cumming et 
al. (2001) suggest that purely radial hydrostatic expansion may 
not be sufficient to explain _rather large Ai/s observed in some 
bursts; see discussion in § 4.3 



It has long been expected that NSs in LMXBs are progeni- 
tors of millisecond radio pulsars (see Bhattacharya 1995 for a 
review). Detection of burst oscillations has certainly bolstered 
the idea that accreting NSs can reach periods similar to those 
seen in millisecond radio pulsars. However, except for one 
source, there is no convincing evidence of millisecond pulsar- 
like - 10 9 G magnetic fields in LMXBs. Only SAX J1808.4- 
365 (Wijnands and van der Klis 1998; Chakrabarty and Morgan 
1998) shows coherent v s = 401 Hz pulsations in persistent emis- 
sion, which is consistent with the fact that a 10 8 — 10 9 G mag- 
netic field channels accretion onto magnetic polar caps and cre- 
ates a permanent brightness asymmetry on the NS (Psaltis and 
Chakrabarty 1999; though, since persistent flux modulation is 
only a few percent, a weaker field could perhaps be sufficient). 
None of the other LMXBs show such pulsations in persistent 
emission, implying that, if they do have magnetic fields, they 
must be weaker than at least 10 9 G. Moreover, the presence of 
frequency drifts during bursts may place even more stringent 
constraints on the magnetic field. If the radial lift-up model de- 
scribed above were correct, it would imply that the ocean and 
atmosphere of the neutron star can make up to fburstA^ ~ 10 
revolutions around its interior. The magnetic field, if present in 
the shearing layer, will be amplified by the shear. A simple esti- 
mate (Cumming and Bildsten 2000) shows that even aw 10 6 G 
magnetic field might be dynamically important in this situation. 
Throughout this work, we neglect the dynamical effects of the 
neutron-star magnetic field, and only briefly discuss possible 
MHD effects at the end. 

In this paper, we abandon the requirement of spherical sym- 
metry, and consider instead the two-dimensional spreading of 
the nuclear fire around the NS surface. We study hydrodynamic 
flows that arise in the burning ocean due to the combination of 
its inhomogeneous lift-up and the rotation of the star. We ana- 
lyze how these flows affect the spreading of fire around the NS 
surface. First we consider the local conditions around the flame 
front at the interface between the hot, burned ashes and cold, 
unburned fuel, exposing the non-trivial effects of rotation and 
viscosity on the propagation of the flame front. Armed with the 
understanding of local conditions, we then construct a global 
scenario for X-ray bursts. 

The burning front during type I X-ray bursts can propagate 
either by deflagration (Fryxell and Woosley 1982b, Hanawa and 
Fujimoto 1984, Bildsten 1995) or by detonation (Fryxell and 
Woosley 1982a, Zingale et. al. 2001). Detonation is possible 
when the timescale for nuclear burning is less than the verti- 
cal sound crossing time. This requires a large column of ac- 
cumulated fuel (~ 100 meters) before the runaway starts, and, 



hence, rather low accretion rates of <J 10 _1i ' 5 Mq yr -1 , i.e., 1 -2 
orders of magnitude less than the observed accretion rates in 
most bursters. However, direct numerical simulations of the 
type performed by Fryxell and Woosley (1982a), and especially 
by Zingale et. al. (2001), if extended into the parameter regime 
relevant for most of the observed sources (i.e., higher accretion 
rate, and, hence, deflagration rather than detonation), are the 
only definitive way to model X-ray bursts. 

In most bursts, the helium/hydrogen ocean burns by spread- 
ing of a deflagration front; its propagation speed is set by the 
rate of heat transport across the front. Fryxell and Woosley 
(1982b) argue that in bursts with quick (< 1 s) rise, which are 
of most interest to us, heat is conducted sideways by convec- 
tion. They give three different phenomenological estimates of 
the speed of the front [Eqs. (3), (4), and the bottom line of the 
left column on page 333 of their paper]. 

1 . The width of the front is equal to the lengthscale of a con- 
vective roll, taken to be approximately the scaleheight h. The 
front speed is then 
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where t„ is the timescale of nuclear burning during the burst. 
This estimate is based on the earlier work of Ruderman 
(1981). 

2. Heat is transported by turbulent convective diffusion, with a 
kinetic diffusion coefficient D = hv c , where v c is the charac- 
teristic convective speed. The front speed is then 
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3. The turbulence scale is much larger than the front width, in 
which case the burning front becomes wrinkled, accelerating 
heat transport. In this case (Williams 1965) 



Vflan 



5 x 10 6 cm s" 



(3) 



Cases 1 and 2, with typical expected t„ ~ O.lsec, h ~ 10 3 cm, 
and v c ~5x 10 6 cm s" 1 , give burst rise times longer than the 
ones observed. Case 3 is, Fryxel and Woosley (1982a) argue, 
the upper limit for the front velocity, but it seems that it is the 
only one of the three estimates which can account for the ob- 
served short rise times of bursts. 

In this paper we present a new mode of heat transport across 
the front which does not depend on a large turbulence scale and 
on wrinkling of the front, and can also explain the observed 
short rise times of type I X-ray bursts. The hot, burned ashes 
have a larger scale height than the cold fuel. This vertical lift- 
up of the ocean as the flame propagates around the star leads 
to a horizontal pressure gradient. Previous work (Strohmayer 
et al 1997, Cumming and Bildsten 2000) only considered the 
situation where the entire ocean is burned and has lifted up, 
and, hence, the horizontal pressure balance is already restored. 
We argue that, as the thermonuclear flame propagates around 
the star, the lift-up of the ocean behind the front drives a dif- 
ferential shear across the front. This shear transports entropy 
from the hot ashes to the fuel, and, hence, propagates the front. 
The speed of the flame depends on the strength of the frictional 
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coupling (e.g., due to convection) between the different layers 
(from top to bottom) of the burning ocean. 
In Section || we estimate analytically that 

(a) In the case of the weak frictional top-bottom coupling, the 
width of the front A is given by the Rossby adjustment radius 
a R : 



A~a R = v / gh hot /f = 2km 



(4) 



2 x 10 14 cm s" 2 / V 10 3 cm 



'hot 



"I 1/2 



2000 rad s" 
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where g is the gravitational acceleration at the NS surface, 
/jhot is the scaleheight of the burned ocean, and / = 257 cos 9 
is the local Coriolis parameter (17 is the angular frequency of 
the NS and ir/2-9 is the latitude). The speed of the front is 



A _ Vg/ihot 

Vflame ^ ~ r, 
tn J*n 

= 2 x 10 6 cm s" 1 



(5) 



2 km 



0.1 s 



which gives a rise time of ~ 0.7 s for burning to spread from 
the equator to the poles [cf. Eq. (g4J)]. 

(b) Let ff r be the timescale of frictional coupling between the 
top and the bottom of the ocean; ff r ~ hh ot /v c if the layers of 
the ocean are dynamically coupled by convection. If the cou- 
pling is strong, i.e., ff r <; ?,„ friction modifies the structure of 
the burning front and its propagation speed. Maximum speed 
is reached when ff r w 1 //, and is given by 
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( ghho 



1/2 



: 3 x 10 7 cm s" 1 
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[Equation ([34|) gives the expression for the burning front 
speed when the frictional coupling strength is arbitrary]. The 
burst rise time is ~ 0.1 s in this case. 

In Section |] we set up and solve numerically a 2-layer 
shallow-water model which contains the essential physics of 
how a stably stratified ocean responds to inhomogeneously ap- 
plied heating. We use this 2-layer model to simulate ignition 
and propagation of the deflagration front, both for weak and 
strong frictional coupling. Our simulations are in agreement 
with the analytical estimates of the front speed in Section |[ 
The shallow-water model also provides us with the theoretical 
estimate of the size of the initial ignition spot on the neutron- 
star surface. However, we are unable to determine the strength 
of friction from the first principles, and we leave it as a free pa- 
rameter of our model. Hopefully, future direct numerical simu- 
lations will address this issue. 

Armed with our understanding of flame propagation, in § |] 
we set out to construct a global scenario for X-ray bursts. We ar- 
gue that, due to rotationally reduced gravity, the thermonuclear 
runaway is likely to begin in the equatorial region. We then 
show that, due to the reduced Coriolis parameter at the equa- 
tor, the flame propagates faster along the equator than away 
from the equator. Thus, a possibly inhomogeneous equato- 
rial belt is ignited at the beginning of the burst rise, and the 



flame proceeds to burn from the equator to the poles. Inside 
the burning front there are strong zonal currents going forward 
relative to the star's rotation, and there are weaker backward 
zonal flows in the cooling wake of the front. A burning inho- 
mogeneity observed as a flux modulation during the burst will 
be trapped in the backward zonal flows, thus the observed mod- 
ulation frequency is smaller than the neutron-star spin. As the 
ocean cools, the zonal flow slows down, and the modulation 
frequency asymptotes to the neutron-star spin frequency. This 
scenario for the burst frequency drift is an extension of the ra- 
dial lift-up scenario proposed by Strohmayer et. al. (1997) as 
modeled in detail by Cumming and Bildsten (2000). 

We speculate that differential zonal currents can be unstable 
to the formation of vortices of the type observed in the atmo- 
spheres of giant planets. These vortices may be the cause of 
oscillations in the X-ray flux in the tail of the burst. In addition, 
if there are vortices trapped in the strong forward zonal cur- 
rents in the burning front, we expect to see a chirp with a large 
frequency span (~ 25 - 500 Hz) during the rise of the burst. Fi- 
nally, we also suggest that the burning front may generate, via 
an MHD dynamo, a magnetic field B ~ 10 9 G, with the ocean 
scaleheight as the characteristic coherence length. This mag- 
netic field can quench the vertical shear in the backward zonal 
flow, and thus may be responsible for the observed coherence 
of burst oscillations in the burst tail. 



2. PROPAGATION OF THE FLAME: ANALYTICAL ARGUMENTS 
2.1. Vertical force balance 

During the burst, the burning material reaches temperatures 
of ~ 2 x 10 9 K, degeneracy in the helium layer is lifted, and the 
burning ocean expands by 10 -40m — many scaleheights of 
the cold pre-burst ocean (Joss 1977, see Cumming and Bildsten 
2000 for the most recent calculations). This vertical expansion 
is very subsonic. However, since only a part of the star is bur- 
ing at a given time, the horizontal pressure imbalance leads to 
non-trivial hydrodynamic flows, as we now discuss. 

Consider a propagating burning front, as illustrated in Fig. [j], 
where the ocean behind the burning front is already hot, while 
the ocean ahead of the front is still cold. Let the ocean, for 
simplicity, reside on a plane, with x and z being the horizontal 
and the vertical coordinates respectively; z is taken to increase 
upward. Since the vertical sound crossing time is much smaller 
than characteristic nuclear burning timescale during the burst, 
?„, vertical hydrostatic equilibrium is always a good approxi- 
mation (we also neglect the vertical component of the Coriolis 
force, see § 2.2). Therefore, 



dp 



= -gp, 



(7) 



where p(x,z) is the pressure, and p(x,z) is the density. From 
Eq. (Q) we see that the separation between constant pressure 
surfaces scales as 1/ p. The hot ocean is less dense than the 
cold one, so the pressure surfaces diverge as they traverse the 
front from the cold into the hotter part, see Fig. [j]. This implies 
that the horizontal pressure gradient dp /dx (indicated by hori- 
zontal arrows in Fig. |lj) increases with height inside the burning 
front. This difference in the horizontal pressure gradient drives 
a shear flow and circulation across the front. 
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Let us quantify the above argument. Let 4>(x,z) be the effec- 
tive gravitational potential per unit mass, dip = gdz. The hori- 
zontal 1 acceleration of a fluid element due to the pressure gra- 
dient is given by 



1 



p(x, (/)) 



dp(x, <j>) 
dx 



Furthermore, 



dp\ 



Therefore, from Eqs. (Kb and (M), we have 



dx 



(8) 



(9) 



(10) 



Eq. ( |To| ) is easy to understand: the height 4>/g of surfaces of 
constant pressure decreases from the hot to the cold part of the 
ocean, so apressm-e is directed from the hot region to the cold re- 
gion (arrows in Fig. [IJ). Differentiating Eq. ( |l0| ) with respect to 
In p at constant x, and using Eq. (0), we get 
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where c s is the speed of sound and 7 is the adiabatic index of 
the gas. As an example, for an ideal chemically homogeneous 

gas 



da-o 



d In p 



0T 

dx 



(12) 



where T is the temperature and p is the mean molecular weight 
of the gas. The term on the right hand side of Eq. ( p"2| ) is the 
forcing term for the vertical shear; it is determined by the hori- 
zontal temperature gradient. 

2.2. Horizontal force balance and flame propagation in the 
ocean without friction 

Let us now write down the horizontal momentum equa- 
tion. In a frame rotating with the star, a fluid element moving 
with velocity v experiences the Coriolis acceleration, a C orioiis = 
— 2f2 x v, where £1 is the angular velocity of the star. In general, 
the Coriolis vector 20 has both vertical and horizontal com- 
ponents in the local frame of the ocean. In what follows we 
shall neglect the horizontal component of the Coriolis vector 
[this is commonly known as the "traditional approximation" in 
geophysics (see, e.g., Pedlosky 1987)]. In doing so, we ne- 
glect two terms in the momentum equation. First, in the equa- 
tion for the vertical force balance (§ 2.1) we neglect the verti- 



cal component of the Coriolis force due to horizontal motion, 
~ flv. As we shall see below, v max ~ (gh) 1 ' 2 , so this compo- 
nent of the Coriolis force is negligible compared to gravity so 
long as SI < (g/h) 1/2 =4x 10 5 (/i/10 3 cm)" l/2 , which is always 
satisfied for neutron stars in LMXBs. Second, we neglect the 
horizontal component of the Coriolis force due to radial mo- 
tion, ~ Clv z , which is small compared to the horizontal Cori- 
olis force due to horizontal motion, ~ flv. This can be seen 



as follows: during the burst, the ocean expands vertically on 
the nuclear timescale t„, and the vertical velocity is v z ~ h/t„. 
Therefore, the horizontal Coriolis acceleration due to radial mo- 
tion is much less then the horizontal Coriolis acceleration due 
to horizontal motion so long as 1 /t„ <^{g/h) x l 2 . This condition 
is always satisfied for type I X-ray bursts. 

With the above approximations, and in the absence of viscos- 
ity (this assumption is alleviated in § 2.3), the horizontal force 
balance is 

dv 



dt 



+/xv, 



(13) 



where v is the horizontal component of the velocity and / = 
2£1 cos 6*z is the vertical Coriolis vector. 

Accreting neutron stars which display type I X-ray bursts are 
typically spinning with frequencies v s of a few hundred hertz. 
The characteristic timescales of the burst light curves are or- 
ders of magnitude greater than the rotation period of the neu- 
tron star: a typical burst rise time is 0.1 - 1 seconds, and a typ- 
ical burst cooling time is 5- 100 seconds. Theoretical calcula- 
tions (see Bildsten 1998 for a review) suggest that the nuclear 
burning timescale during a typical burst is t„ ^0.1 sec, still 
much greater than \/v s . Therefore, during the burst the ocean 
flow must be in quasigeostrophic equilibrium (Pedlosky 1987) 
everywhere except for the immediate vicinity of the equa- 
tor. Quasigeostrophic equilibrium implies that inertial external 
forces acting on a fluid element in the horizontal direction — in 
particular, due to the horizontal pressure gradient — are almost 
exactly balanced by the horizontal component of the Coriolis 
force. 

One can see this by noting that the second term on the right- 
hand side of Eq. (|l3]) is dominant. We can rewrite this equation 

as 

1 ' ^ 'r/xf. (.4) 



Af X 3 
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dt 



To zeroth order, 
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-fx d p 



(15) 



Substituting this expression into the right-hand side of Eq. ( |14| ) 
we get the first-order correction to the velocity: 



V — Ageostrophic 
1 



+ V. 



- f X £Z press 



ageostrophic — 

1 da 



pressure 



f- 



dt 



(16) 



The first term in Eq. (JT^) represents the dominant geostrophic 
flow in the direction along the front line and perpendicular to 
the pressure gradient. The second term is the ageostrophic com- 
ponent of the velocity, which is perpendicular to the front line. 
Differentiating Eq. ( |l6j ) with respect to In p and using Eq. (|LT 
we get 
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and 



dv. 
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1 d da-a 



d In p f 2 din p dt 



(17) 



(18) 



Equation ( |17[ ) is commonly referred to as the thermal wind 
relation in the meteorology/geophysics literature. It implies that 



1 A horizontal surface is a surface of constant <j>; in plane-parallel geometry with constant g, setting </>=const is equivalent to z=const. 
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the difference in geostrophic velocity <5v ge0 strophic between the 
top and the bottom of the ocean somewhere in the middle of the 
front is related to the difference in p/p between the hot and the 
cold parts of the ocean: 



Sv, 



geostrophic 



top _ bottom 

geostrophic v geostrophic 



(P/P)h0t-(P/P) 



cold 



gh 



hot 



/A 



/A' 



(19) 



where /ihot is the scaleheight of the ocean behind the front which 
has just undergone a thermonuclear runaway, A is the width of 
the front, and we used the fact that /ihot ^coid- 

Likewise, we can use Eq. ( |l8| ) to estimate the characteristic 
ageostrophic shear <5v ageostrophic = v a ° g p eostropmc - v**S£oi«c some- 
where in the middle of the front. Since the front passes through 
a given fluid element on the timescale of order t n , we replace 
d/dt by 1 jt n in our order-of-magnitude estimates. Then 



dv, 
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■ Sv 



1 da B 
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gh 



hot 



f-t„ dlnp /2f„A 



n. (20) 



Here n is the unit vector perpendicular to the front line. 

The quantity <5v a geostrophic is tne characteristic speed of the 
shear flow perpendicular to the burning front. How does it re- 
late to the flame propagation speed? Generally, flames propa- 
gate by transporting entropy from the hot, burned material to 
the cold, unburned fuel, and their speed is set by the time it 
takes to transport this heat across the width of the front. In 
the case of slow, laminar deflagration fronts (see, e.g., Bildsten 
1995) heat is transported by conduction, while for convective 
fronts discussed by Fryxell and Woosley (1982b), heat is ad- 
vected by the turbulent motion of the fluid and then mixed into 
the fuel. In our case, however, it is the shear flow, £>v a geostrophic, 
that moves the hot material ahead of the front (this happens at 
the top of the burning ocean), and pulls the cold fresh fuel into 
the burning front (this happens at the bottom of the ocean). 

We assume that vertical mixing occurs within the front, and 
that this mixing transports heat between the hot burning fluid 
moving forward on the top and the cold fluid on the bottom. 
This mixing could be either due to convection, or due to var- 
ious shear instabilities. Indeed, one-dimensional numerical 
simulations, starting with Joss (1977), show that X-ray bursts 
are strongly convective, with convective overturn timescale of 
~ 10~ 3 s. If convection is present within the front, it will ef- 
ficiently mix entropy in the vertical direction. However, while 
in one dimension the hot burned fluid has no place to go but 
to mix with the cold unburned fuel, in more than one dimen- 
sion convection may be quenched by the sideways propagation 
of the hot material on top of the cold material 2 . In addition to 
the possibility of convection, the geostrophic flow within the 
front possesses strong vertical shear, with velocity difference 
of Sv ~ Vg^hot across a scaleheight [cf. Eq. ( |l7| ) and (|25|)], 
which corresponds to the Richardson number of order unity. 
The front with such strong shear may be unstable to the Kelvin- 
Helmholtz instability. Finally, the interface between the shear- 
ing layers is not exactly horizontal, but has a slope of order 
fthot/A <; 10~ 2 . Such flow may be unstable to the baroclinic 
instability (Fujimoto 1988, Cumming and Bildsten 2000). The 
nonlinear development of these or other instabilities might re- 
sult in efficient vertical mixing. The presence of efficient ver- 
tical mixing within the front is the most uncertain part of our 

2 We thank Lars Bildsten and Andrew Cumming for making this point. 



model; it must be addressed directly by future numerical simu- 
lations. However, given the large amount of shear in the burning 
front, such mixing is not at all unreasonable. 

We therefore believe that the vertical thermal mixing 
timescale f m j x i ng within our geostrophic front is much smaller 
than t„ (e.g., the convective overturn timescale is 10~ 3 sec). In 
the foregoing, we assume that the entropy advected forward by 
the ageostrophic flow is quickly mixed in the vertical direction. 
If the thermal mixing timescale f m j x i ng within the front is non- 
negligible, then the general argument laid out in this paper is 
unchanged, except that all estimates should use t„ + f m i x i ng in- 
stead of just t„. 

If entropy mixing is efficient, then the front propagation 
speed is the velocity with which entropy is transported, i.e., the 
characteristic speed of the ageostrophic shear flow: 



Vflame ~ Sv, 



gh 



hot 



ageostrophic ^ ro A ' 



(21) 



The width of the propagating front is (e.g., Fryxell and Woosley 
1982b) 

A ~ Vflamefn- (22) 

Substituting this into Eq. (|T]), and solving for the front speed 
and width, we get 
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(23) 
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(24) 



for g = 2 x 10 14 cm s~ 2 . From Eq. ( p^ ) the characteristic 
geostrophic velocity inside the front is of order the gravity wave 
speed in the hot material: 



Ageostrophic ~ (,?^hot) 



1/2 



(25) 



Note that the width of the front is equal, to the order of magni- 
tude, to the Rossby adjustment radius (Pedlosky 1987). 

2.3. Horizontal force balance and flame propagation in the 
ocean with friction 

In the previous subsection we have assumed that there are 
no viscous forces acting on fluid elements, and hence that the 
top and the bottom of the burning ocean are free to slip past 
each other. However, there is likely to be some vertical mix- 
ing within the front, due to turbulence driven by convection or 
by strong shear (see the discussion in the previous subsection). 
This mixing will exchange momentum between the top and the 
bottom of the ocean; we must, therefore, consider what effect 
viscous-type friction within the ocean will have on the propa- 
gation of the burning front. 

Let us return to the horizontal force balance equation ( |l3| ) 
and include viscosity: 



dv 
dt 



+/x v- 



(26) 
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where we have added to the right-hand side the acceleration 
of the fluid element due to viscous stress. We now consider 
the difference in pressure between the bottom and the top of the 
burning ocean, where we take the top to be roughly one scale- 
height above the bottom. From the previous equation we have 



^op ^ottom 

pressure pressure 



dv 



top 



dv 



bottom 



dt 



+ /x(v 



top 



dt 

" v bottorrJ v fl visc( 



_ ^-bottom \ 
a viscous'' 



(27) 



We assume that effective viscosity dynamically couples the top 
and the bottom of the burning ocean, and model this coupling 
by introducing a linear friction term: 



/^"P _ / jbottom f-ff _,7 "v 

fl viscous "viscous — . v v top ^bottom.)- 
ffr 



(28) 



If friction is mediated by turbulent eddies, then 1 /ft = (3\ v tur b /h, 
where v tur b is the characteristic velocity of the largest turbulent 
eddy, and (3\ is a number of order unity. On the left-hand side 
of Eq. (Evh we use Eq. (ph to get 



-top -bottom 

pressure pressure 



d\np 



d( P /p) 



dx 



h 



hot . 



(29) 



Finally, as in § |2.2| , we replace in ( |27| ) the time derivative d/dt 
by rj/t„, where r\ is a number of order unity. With the above 
simplifications, the shear flow i5v = v top - bottom obeys 



/iho 



f x Sv+ ( — + — ) ('li- 
fe t n 



Scalar multiplication of the above equation by Sv yields 

^hot . 



-{dv ■ h) ■ 



(30) 



(31) 



while application of Pythagorean theorem to Eq. (|30|) yields: 

2 



(Sv) 2 



gh*, 



l 



A J / 2 + (l/fe+r?A„) 2 " 



(32) 



Therefore, the component of the shear flow across the front is 
just 

,^ „ ghhot(l/tfr+r)/tn) 



[/ 2 + (lA fr + 77/f„) 2 ]A- 



(33) 



As argued in § 2.2, the speed of the front is determined by the 
speed of the shear flow across the front line, i.e., Vfl ame ~ Sv ■ h. 
Then, using A ~ Vfiamef« in Eq. ([33]), we get 



Vfla 



gh 



hot 



l/fe + 7?A„ 



t„ f 2 + (\/t b + r)/t n ) 2 



1/2 



(34) 



The typical parameter values relevant for the burning ocean 
are g = 2 x 10 14 cm s" 2 , hhot = 10 3 cm, t n = 0.1 s, and / = 
2000 rads" 1 . For purposes of illustration, we take the eddy 
velocity to be equal to the typical convective velocity infered 
from 1-D models of x-ray bursts v tur b = v c = 5 x 10 6 cm s" 1 , and 
set the frictional timescale ff r to /fhot/vmrb- We then obtain the 
numerical value of v name ~ 200 km s -1 . This corresponds to a 
burst rise time (equator to pole spreading) of ~ 0. 1 s. 



In Figure || we show the speed of the flame front Vflame as 
a function of the strength of the dimensionless frictional cou- 
pling constant, (t/ r f) ■ The dots represent the results of our 
model simulations, which are discussed in the next section. 
The continuous line is a fit made using the functional form of 
Eq. (Q) and adjusting the values of and t„ to match the 
simulation. The agreement between the model simulations and 
the analytical formula for the speed propagation is very good, 
which means that the analytical estimates capture the physics 
contained in the problem quite well. 

In Figure ||, three regimes of the front speed are evident. 
When frictional coupling is very weak (i.e., ff r ^> t„ and ff r ^> 
1 If) the front speed asymoptotes to a constant value given by 
Eq. (E4), i.e., we recover the result of § 2.2 for frictionless flame 



propagation. When friction is important (i.e., ff r < ?„), there are 
two distinct regimes, depending on whether the coupling fre- 
quency 1 /ff r is smaller or greater than the Coriolis parameter /. 
When 1 /ff r -c / (friction is small), then 



Vila 



g^hot l/tb + Tj/tn 



ftn 



f 



1/2 



while if friction is strong (1 /f tr ^> f) 



Vflair 



ghhottfr 



1/2 



(35) 



(36) 



For a fixed frictional coupling timescale ff r , the front speed is 
always a decreasing function of /; for a fixed / the front speed 
attains a maximum value when friction is acting on the rotation 
timescale, 1 /ff r = /: 



max 
flame 



g/?hot 
ftn 



1/2 



(37) 



Comparison with Eq. (|24|) shows that the maximum speed 
is a factor of {ft n ) l l 2 ~ 15 larger than the front speed without 
friction. The increase of the front speed when friction is intro- 
duced into the ocean can be qualitatively understood as follows. 
Without friction, the dominant component of the fluid velocity, 
Vgeostrophic, is exactly parallel to the front, and the cross-front 
velocity, (5v ageostrop hic, is a factor of ft„ smaller. The presence 
of friction modifies geostrophic balance so that the geostrophic 
shear flow has a component perpendicular to the front. How- 
ever, if friction becomes too large, it will suppress the shear, and 
the front will stall. The cross-front component of geostrophic 
shear is maximized when friction and Coriolis force have the 
same magnitude, 1/fe w /. 

In the above analytical estimates we neglected the effect of 
the magnetic field that may be present on the neutron star sur- 
face. What is the magnitude of a B field that would have dy- 
namical consquences? This B field would have to alter the 
leading-order geostrophic balance [cf. Eqs ( |l3| ) and ([26|)], i.e., 
^dynamical ~ (47rpv 2 eostrophic ) 1 ' /2 . If viscous coupling is negligi- 
ble, Vgeostrophic ~ (?^hot) 1/2 , and 



B, 



dynamical ' 



10 12 G 



p 



10 6 g cm -3 



1/2 



10 3 cm 



1/2 



(38) 



i.e., magnetic energy density needs to be comparable to the 
hydrostatic pressure in order to affect geostrophic flow. Con- 
versely, if convective viscosity is such that Vflame = v Sl & < then 
Vgeostrophic ~ gh/fA ~ {gK ol / ft n Y I 2 , and the above estimate is 
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a factor of (/f„) 1//2 ~ 15 smaller. Even in order to affect the 
first-order correction to the geostrophic balance, i.e., the cross- 
front circulation, the magnetic pressure has to be comparable to 
/ 3V ageostrophic/2- The smallest ageostrophic velocity is attained in 
the case of no friction, v ageostrophic ~ (gh hol ) l/2 / ft n [cf. Eqs @ 
and (p4[)], which still requires Bdynamicai <; 10 9 G in order to af- 
fect the details of the cross-front circulation. As discussed in 
the introduction, it is unlikely that accreting neutron stars in 
LMXBs possess surface magnetic fields of this magnitude. 

The argument above applies to the instanteneous effect of 
magnetic field on the flow, and does not take into account the 
nonlinear evolution of magnetic field in the (possibly turbulent) 
shear flow, and the back-reaction of the evolved field on the flow 
(see, e.g., Kluzniak and Ruderman 1998, Spruit 1999, Cum- 
ming and Bildsten 2000). Detailed consideration of this issue 
is beyond the scope of this paper. In Section 4 we only briefly 
speculate on magnetic field generation by an MHD dynamo in 
the burning front. 

In the next section we construct an explicit 2-layer shallow- 
water model for flame propagation, which we designed to in- 
clude the essential physics described in this section. Our nu- 
merical simulations confirm the analytical results presented 
here. 

3. PROPAGATION OF THE FLAME: NUMERICAL MODEL 
3.1. Motivation and basic description of the two-layer model 

In this section we outline one particular approach to mod- 
eling hydrodynamic flows in a stratified plane-parallel atmo- 
sphere 3 in the presence of localized heat sources and sinks 
(thermal forcing). The thermodynamic state of a parcel of fluid 
in the atmosphere will be described using its temperature, T, 
pressure, p, and potential temperature, 8 = T(po/p) K , where 
po is some constant reference pressure, and n is related to the 
adiabatic index of the gas, 7, as k= 1- I/7. Potential tem- 
perature can be interpreted as the temperature which a parcel 
would have if it were adiabatically brought from pressure p to 
the reference pressure po (Holton, 1992). It equals to exp(s), 
where s is the specific entropy of the parcel. In a stably strat- 
ified atmosphere potential temperature should increase or be 
constant with height. Another property which follows from 
the definition of 9 is that for fixed pressure it is proportional to 
the specific volume of the gas: 9 cx ^/? (I_K) . For the particular 
case k = 1 which will be utilized below, potential temperature 
is equivalent to the specific volume of a parcel of fluid. 

Now we can begin to evaluate the qualitative effects of heat- 
ing on an atmospheric column. From the heat equation we have: 



d\n9 J 



dt 



c p T ' 



(39) 



where J is the net rate of heating (or cooling) per unit mass, 
which, in our case, is the difference between the nuclear energy 
generation rate and the radiative cooling rate. The effect of 
heat input, therefore, is to change the potential temperature of 
the element of fluid. Such a change, coupled with maintenance 
of vertical hydrostatic balance in the column leads to hydrody- 
namic circulation. As a simple example of how such a circula- 
tion arises consider two columns of gas of equal column density 
but of different potential temperature (9\ > 62) positioned next 
to each other (see Fig. |3|a). Each of the columns is assumed to 
be in vertical hydrostatic equilibrium. This situation can arise 

3 Henceforth we use the terms "atmosphere" and "ocean" interchangeably. 



if column 1 is strongly heated, so that its entropy is instantly in- 
creased relative to the entropy of column 2. This configuration, 
however, is not stable: the fluid with higher entropy will end up 
on top of the fluid with lower entropy, and the atmosphere will 
become stably stratified (see Fig. ||b). The dynamics of this ad- 
justment can be qualitatively described as follows. Initially, let 
the tops of the columns be at the same outside pressure; the bot- 
tom pressure then will also be the same since the column densi- 
ties are equal for both columns. At higher altitude the pressure 
in column 1 is larger than in column 2 at the same height be- 
cause the density in column 1 is lower than that of column 2 
(cf. Fig. [j]). At first the fluid will flow from column 1 towards 
column 2 along the top, and not along the bottom where the 
pressures are the same. As the high-entropy fluid accumulates 
on top of column 2 the pressures at the bottom no longer match, 
and there will be a flow of fluid from column 2 towards column 
1 along the bottom as well. This circulation will continue un- 
til the equilibrium state shown on Fig. |]b is achieved. Note 
that we assumed this equilibration is happening adiabatically, 
so that the fluids of different entropy maintain their identity and 
in the end do not mix but stratify. 

The entropy stratification described above suggests a useful 
analogy for the effect of a heat source in the atmosphere. Since 
higher entropy fluid ends up on top of the lower entropy fluid, 
we can treat entropy (or potential temperature) as a general- 
ized vertical coordinate. The action of the heat source is then 
to pump fluid from a low entropy level to a high entropy level 
while conserving mass. The isentropic levels generally do not 
coincide with the levels of constant height, and there will be a 
circulation of fluid. If the motion outside the heat source is adi- 
abatic, it is impossible for the fluid to cross the isentropic level 
again, and the motion will be along the isentropes. 

Henceforth, we will adopt a particularly simple, yet powerful 
model which uses the equation of state with n = 1 (7 — > 00). The 
fluid with such an equation of state is quasi-incompressible: it 
does not allow adiabatic motions which involve compression 
and rarefaction, and has an infinite speed of sound. However, 
for k = 1, the potential temperature is equal to the specific vol- 
ume of the fluid, and heating will lower the density of the fluid 
while conserving its mass. 

Following Ooyama (1969) we consider an atmosphere con- 
sisting of two layers of incompressible fluid of different density, 
with density ratio e = pa/pi < T The lighter fluid is assumed 
to be in the top layer (layer 2). A fluid element is allowed to be 
in one of the two possible density states, p\ or p 2 - The heating 
is assumed to be located at the interface between the two lay- 
ers (which is not necessarily a surface of constant height); the 
physical effect of the heating is to convert the "cold" dense ma- 
terial of density p\ into the "hot" light material of density P2, 
transferring the fluid from the lower (cold) to the upper (hot) 
layer. The equations of motion for the two layers are then given 
by the shallow-water system of equations (Pedlosky 1987) in- 
dividually formulated for each layer. The crucial distinction 
of our model is the interlayer coupling source/sink term in the 
mass continuity equation, and the friction between the layers: 



dh. 

_i + V| r /j 1 vi=-e(e + -G_), 

^+v r /i 2 v 2 = e + -e_, 

dv\ x 

+ vi • V| 1 vi x = -gV x (h 1 + eh 2 ) + fv\ y 



(40) 
(41) 
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dt 



--(Q- + Pf)(v lx -v 2x ), (42) 
+ T 2 • V 1 1 V2x = -gV x (h i+h 2 ) + fv Zy 

- ^P^(v2,-v lx ), (43) 



• V||Vi, = -fVU-j^ (fi- + A*/) (yiy-V 2 y), (44) 
+ V 2 • V| | V 2> . = "/V2x - ^7^ (V2 - V " Vll) ' (45) 

cf. Eqs (3.5)-(3.7)of Ooyama (1969). Here, h\2 are the heights 
of the cold and hot layers respectively. We consider the dynam- 
ics on an x,y plane with V 1 1 = (d/dx,d/dy). To keep the prob- 
lem tractable we will allow variation only in x direction, so that 
we have a slab geometry in x with possible tangential velocity 
v y . Depending on the case under consideration, the plane of the 
simulation will be oriented in various ways on the neutron star 
surface. The heating rate Q = J/g represents the flux of ma- 
terial from the lower to the upper layer, and has dimensions of 
velocity (for the derivation of Eqs. ([MJ — see the appendix). 
The quantities Q+ and Q- are defined as follows: Q+ = Q, Q- = 
if the heating rate Q > 0, and g_ = -Q, Q + = if Q < 0. The 
terms on the right-hand side of the above equations which are 
proportional to p f(v\ - V2) and Q±(y\ - v 2 ) represent the friction 
between the two layers, the former due to (perhaps convective) 
viscosity, and the latter due to momentum-conserving transport 
of material from one layer to another. 

In general, the heating rate depends on the temperature at the 
point where the heat source is applied. Ooyama (1969) argues 
that for the shallow water model it is reasonable to assume that 
the temperature of the fluid can be represented by its height. 
His argument is roughly as follows: an adiabatic temperature 
rise in a normal fluid happens when the flow is converging; but 
a converging flow in a shallow-water fluid results in a rise of 
its height. Likewise, an increase of the temperature due to heat 
input corresponds to addition of fluid to a shallow-water layer, 
and hence corresponds to an increase of its height. Thus the 
temperature increase in a normal fluid, both due to compression 
and due to the heat input, corresponds to the height increase in 
a shallow-water fluid. In our model, we shall assume that the 
temperature at the interface where the heating occurs is set by 
the height of the second, light layer through a simple relation- 
ship T = ghijcp. 

In the following three subsections we consider how our 
model responds to localized sources of heat. This will help us 
develop intuition which will be useful in analyzing the ignition 
and spreading of the burning front. 

3.2. Localized heating in a nonrotating atmosphere without 

friction 

In a nonrotating atmosphere the effective Coriolis parame- 
ter / is zero, and this may be interpreted as the burning in 
a region along the rotational equator of a neutron star. Con- 
sider a part of the atmosphere of a large lateral extent with the 
two-layer fluid described above in the equilibrium configura- 
tion with, for example, h\=hi initially. The fluid is allowed 
to freely leave and enter the domain through the sides. For 
a delta-function heat source Q = QqS(x) and zero friction be- 
tween the layers ([if = 0) we can find a steady state solution 
that will be achieved after an initial transient. The solution 
will have a uniform divergent flow at the upper layer, satisfying 



hiVix = 2o(@( x ) _ 1 /2), and a uniform convergent flow at the 
bottom layer, with h\v\ x \ = -eQo{<d{x)- 1 /2), where 0(x) is a 
step-function. The equilibrium value of /12, after the transient 
flow traverses the domain of interest, is larger by ~ Qo/\/gh 
than the initial thickness of the upper layer, while the layer 1 is 
contracted by similar amount. The main property of this solu- 
tion is that, even though the heat source is continually operat- 
ing, the temperature (or thickness) of the layers does not change 
in time after the initial transient has died out. The heating drives 
opposing winds in both layers, and the heating energy goes into 
raising the fluid from the bottom to the top. For small heating 
rates, Qo/h <; \fgh, the shearing flow between the two layers is 
stable to Kelvin-Helmholtz instability due to buoyant stratifica- 
tion. 

This solution is valid until the outflow has reached the end 
of the lateral extent of the atmosphere, at which point layer 2 
will grow at the expense of depleting layer 1. Such lowering 
of layer 2, i.e., the presence of high-entropy fluid at a lower 
height, represents simple heating of the whole atmosphere that 
we eventually expect from a heat source. Therefore, in order 
to have the heat source continuously increase the temperature 
locally, as is required for a local thermonuclear runaway, one 
needs a way of containing the outflow over some spatial length 
scale. In the absence of lateral boundaries, this feat is accom- 
plished by rotation of the star, friction, or both. 

3.3. Localized heating in a rotating atmosphere 

General considerations. The main new feature introduced 
into the dynamics in a rotating atmosphere is that when friction 
is absent, all disturbances do not relax to a state of minimum 
potential energy, in contrast to a nonrotating atmosphere. In- 
stead, disturbances tend to relax rapidly (in a time of the order 
of one rotation period) to a state of geostrophic balance, where 
all pressure gradients are balanced by the Coriolis force acting 
on fluid moving along the isobars. A famous example is known 
as the Rossby adjustment problem (Gill 1982): a uniformly ro- 
tating fluid of constant density is released with an initial step 
in its height distribution. A transient ensues, in which the fluid 
starts to spread. The spread is slowed down, however, by the 
action of the Coriolis force, and the fluid oscillates around the 
equilibrium configuration, radiating gravity waves. The end re- 
sult is a state of equilibrium with nonzero kinetic energy in the 
flow such that the adjusted pressure gradient is balanced by the 
Coriolis force. A simple way of understanding the spatial scale 
on which the equilibration occurs is to consider the momentum 
equations for a single incompressible layer of constant depth 
Ho with a perturbation h = Hq + Sh: 



dv, 

dt 

dv y 

dt 



d 

-g-7^(H + 6h) + fv y , 



—77 = -fv x . 



(46) 
(47) 



Equation (^) can be integrated to give v y = fAx, i.e., the 
transverse velocity of the fluid element is proportional to the 
displacement from its initial position. Substituting this into 
( fig ) with the assumption of geostrophic balance, and esti- 
mating the pressure gradient as gHo/Ax, we obtain the char- 
acteristic length over which the pressure gradient is spread: 
Ax = \/gHo I f. This is the Rossby radius of deformation for 
a rotating atmosphere of scaleheight Hq. After the transient, the 
initial step in height of the fluid will transform into a gradual 
slope over the lengthscale of Rossby radius. 
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The above argument is valid only if the horizontal displace- 
ment of the fluid element, Ax, is of the same order of magnitude 
as the lengthscale over which the pressure gradient is spread. 
This is true only for large initial perturbation, 8h ~ Hq. How- 
ever, the conclusion that the pressure gradient is spread over 
the Rossby radius turns out to be valid even if the initial per- 
turbation is not large. For a small perturbation, the spread of 
the pressure gradient does not involve actual motion of a fluid 
element over Ax, but is communicated via gravity waves. In 
the presence of rotation, the gravity waves are dispersive, with 
dispersion relation uj 2 = f 2 + c 2 k 2 , where c g = y/ 



o. For wave- 
lengths shorter than the Rossby radius, gravity waves are unaf- 
fected by rotation and leave the region of disturbance, while for 
wavelengths longer than c„/f the propagation speeds become 
very small. Only the long-wavelength component of the initial 
perturbation is left behind, and even for small initial perturba- 
tions the geostrophic ballance is eventually established on the 
lengthscale of the Rossby radius, ^gHo //. We refer the reader 
to Gill (1982) for further details. 

Localized heating in a rotating atmosphere without fric- 
tion. We now study the flows in a rotating atmosphere forced 
by a localized heat input. We concentrate on the top layer of 
the isentropic model discussed above. If we assume that the 
density contrast between the two layers is large enough (e <C 1) 
then the height of level 1 does not change appreciably during 
the evolution. In that case, the linearized equations for the top 
layer of thickness h = Ho + Sh become: 



9h d , 
-^+H —v x = Q(x,t), 

9v x 
dt 
9v y 



d 



(48) 
(49) 
(50) 



We will further refer to models tracking only the top fluid layer 
as one-layer models. We find that they contain most of the es- 
sential physics of the two-layer system. If the timescale of heat- 
ing is much slower than the time it takes gravity waves to prop- 
agate over the characteristic Rossby radius (i.e., 1 //), at each 
moment in time the system will be in approximate geostrophic 
balance. We thus omit the time derivative in ([49|); this also au- 
tomatically filters out transient gravity waves. The system of 
equations can now be combined to yield a single equation for 
the height of the layer: 



-) d h , , 
R dx 2= ' Vf • 



(51) 



This is a forced Klein-Gordon equation, where the source 
term is provided by the total heating at a particular point. 
The Green's function for the left hand side of the equation 
is G(jc|0 = 1 /2aRe~\ x ~'>\/ aR . Therefore, for a constant delta- 
function heating Q(x,t) = Qo6(x) the fluid follows a growing 
vortex sheet solution: 



Kx,t) = -z — e 
2a R 



Go 



v^,?) = signW^e-W/* ! 
v y (x,t) = -sign(x)^fe-^ a «. 



(52) 
(53) 
(54) 



A few points should be made about this solution. It shows that, 
for constant heating, the height of the layer (or its temperature) 
grows linearly in time in a region of characteristic size equal 
to the Rossby radius around the heat source. Aside from non- 
linear effects, the size of the area affected by heating does not 
change as the height increases. Therefore, the pressure gra- 
dient increases linearly with time. The velocity transverse to 
the pressure gradient (geostrophic velocity v v ) also grows lin- 
early with time. If we were to solve the same problem with 
cylindrically symmetric localized heating, the geostrophic flow 
would create a vortex-like circulation around the source of heat- 
ing. When Qo > the heated region is the area of high pres- 
sure, and hence the sense of geostrophic velocity in the layer 
is anticyclonic (opposite to the sense of rotation of the star), 
and cyclonic in a low pressure region for the case of cooling. 
The discontinuity in velocity at x = is an artifact of the delta- 
function forcing. Since realistic forcing is distributed over a 
finite area, the velocity should smoothly go to zero at the center 
of heating. Another feature of the solution is the presence of 
a nonzero flow away from the source of heating (and the op- 
posite for cooling). This flow is not due to pressure gradient 
as such, but rather due to a time rate of change in the pressure 
gradient as the central pressure rises or falls. Such a flow is 
known as the isallobaric wind (Holton, 1992) and is not itself 
in geostrophic balance. However, this wind plays an important 
role in the adjustment to balance as can be seen from Eq. (^0|). 
In order for the geostrophic velocity to change in response to 
modified conditions, there must be a nonzero ageostrophic ve- 
locity; cf. Eq. (18). Since all quantities in the above solution 



exponentially decay at distances larger than the Rossby radius, 
the effects of heating are thus localized due to rotation. 

Localized heating in a rotating atmosphere with strong fric- 
tion. When there is (turbulent) viscosity in the system, fric- 
tional forces need to be added to the right-hand side of Eqs. 
(pb and (Bob: 



dh d „, , 
-^+H Q —v x = Q(x,t), 

dv. x 
dt 

8Vy Vy 

-7- =-fv,-— ■ 

at f fl - 



= - 8 dx' 1 + fV > 



(55) 
(56) 
(57) 



Here l/ff r is the coefficient of frictional drag (cf. § 2.3). For 
l/ f fn/ l/t„ we can neglect the time derivatives on the left- 
hand side of Eqs. (|56|) and (|57|). We then get 



dh d 2 h 
- + D— 2 =Q { x,t), 

_D dh 

x Ho dx ' 

Vy = -ft&Vx, 



where 



D 



t&gHo 



(58) 

(59) 
(60) 

(61) 



l+(/ffr) 2 ' 

Equation (|8j) is a diffusion equation with the source term 
Q and the diffusion coefficient D. For a localized heating 
Q(x, t) = Q()S(x) switched on at t = 0, the solution is given by 



h(x,t) = Q 



4ttD 



F(x/2VDt) 



(62) 
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where 



F(q) = q 



2' 



l dqi 



(63) 



of a star with 250 Hz rotation frequency at 45° latitude. First 
we discuss simulations for the case without interlayer friction 
[M/ = 0in(@-||)]. 

The initial height perturbation introduces potential energy 
about 2/3 of which is radiated into gravity waves, and the rest 
ends up in a geostrophically balanced flow which is established 
within a rotation period and has a characteristic lengthscale of 
a Rossby radius. As expected, the initial stage of growth is sim- 
ilar to the delta-function response described above: the vortex 
intensifies while maintaining its width. Since in the beginning 
of the runaway the heating is a strongly increasing function 
of temperature, the dominant contribution to vortex intensifi- 
cation comes from the fastest-growing central part. The heat- 
ing function (Q) goes through a peak during the runaway at 
about 9 x 10 s K for our parameters before cooling catches up 
with heating at about 2 x 10 9 K. When the central temperature 
moves past the peak heating, and the tails of the vortex start un- 
dergoing the fast part of the runaway, a double-peaked structure 
of the heating function appears as seen in fig. [|. As the heat- 
ing peaks separate, the vortex begins to spread and the structure 
of a burning front appears: a localized area with the runaway 
heating which moves into the unburned material at a constant 
speed. The burning is accompanied by a rise in layer thickness, 
so that there is a substantial pressure gradient in the direction 
of motion of the front. In the absence of friction, this pressure 
gradient is closely balanced by the Coriolis force acting on the 
flow generated parallel to the front. As discussed in Section |[ 
the unbalanced ageostrophic component is present and neces- 
sary to drive the change in the much larger geostrophic flow: 
dv g eostrophicM = -/v ag eostrophic • It is the ageostrophic flow (typi- 
cally smaller in magnitude by a factor of order 1 / ft„ ~2x 10 2 
compared to its geostrophic counterpart) which is responsible 
for the front propagation. 

3.5. Structure and speed of the front 

The temperature at a location ahead of the front increases 
due to two effects: influx of the ageostrophic wind of hot mate- 
rial across the front, and local thermonuclear energy generation. 
The local heating rate due to thermonuclear reactions, which 
sensitively depends on the temperature, will start to grow if the 
temperature perturbation is sufficient to push the fluid parcel 
into the runaway regime. We are assuming that the temperature 
and column depth in the fluid are such that by raising the tem- 
perature sufficiently the runaway will start. In this case, since 
the temperature well inside the burning front is high enough 
for a runaway, we are guaranteed that even if the fluid before 
the front is colder than the point of marginal stability, it will 
cross into the runaway as the front approaches. When the fric- 
tion is absent, the fluid motion across the front is driven by the 
time rate of change in the pressure gradient: v x = -(g/ ' f 2 )^^h 
[cf. eq. (PQ)]. As the front approaches, the pressure gradi- 
ent increases in magnitude (it is negative for a front moving to 
the right) and reaches the maximum at the point of peak heat- 
ing. The fluid is thus pushed through the front, igniting material 
ahead, and the pressure gradient is performing work to acceler- 
ate this fluid to the geostrophic velocity. The width of the front 

4 The heating function that we use is representative of the conditions in a type I burst, but is by no means complete. In particular, we neglected the electron 
screening of the 3a reaction (Fushiki and Lamb, 1987), as well as further energy release due to nuclear evolution beyond carbon. We find that the physics of front 
propagation is independent of the particular heating function as long as this function has a rise, a single peak, and decays at high temperatures. 

s For our numerical method we use second order accurate finite-differencing of model equations. This allows us to cut down on numerical diffusion that can 
artificially increase front speeds. Since the underlying equations support fast moving dispersive gravity waves, special attention is paid to boundary conditions. To 
allow for integration times longer than the gravity wave crossing time, we introduce an absorbing boundary layer (Romate 1992) that dissipates the gravity waves 
of wavelength smaller than the size of the boundary layer. To prevent reflection of longer wavelengths we add the non-reflecting condition that eliminates backward 
going characteristics at the boundary (Thompson 1987). 



The height (i.e. temperature) and velocity perturbations are con- 
centrated within the diffusion length Sx = \[Dt. 

Implications for burst ignition lengthscale. What can we 
learn from the atmospheric response to localized heating? We 
saw that for a frictionless rotating atmosphere, the temperature 
perturbation is confined to the Rossby adjustment radius, hence 
the ignition of the burst is likely to happen on that scale. When 
there is strong friction, the temperature perturbation is confined 
to the diffusion lengthscale. In this case, we therefore expect 
the ignition to happen on the scale <5jq g nition ~ \/Dt n , where D 
is the diffusion coefficient given by Eq. (p5l|), and t n is the char- 
acteristic nuclear burning timescale. As we saw in the previous 
section, and as will be confirmed by simulations in the next sec- 
tion, these ignition lengthscales are of the same order as front 
widths in respective cases; cf. Eqs ( f22[ ) and ([34]). The only dif- 
ference is that in the case of the front propagation the important 
scaleheight is that of the hot part of the ocean, whereas for the 
ignition the scaleheight Ho is that of the unheated, cold ocean. 

3.4. Burning front propagation 

When we numerically solve the system of equations (|40|)- 
( p5[ ) with a temperature-sensitive heating function, we find that 
a local runaway can develop into a propagating burning front 
solution. In order to simulate conditions relevant to the case of 
a neutron star atmosphere during a burst we consider heating 
due to 3a helium-burning reaction and one-zone cooling due to 
black-body radiation (see, e.g., Cumming and Bildsten 2000) 4 : 



Q = 5.3 x 10 21 erg g" 1 s" 1 ^ exp f — 



T 3 

i 8 



-44 \ 
7W 



acT 
3ny 2 



(64) 



Here, p$ is the density in units of 10 5 g cm" 3 (which we evaluate 
including degenerate corrections), Y is the helium abundance 
obeying dY jdt = 63 Q /(5.84 x 10 17 erg/s), 7g is the temperature 
in units of 10 8 K, and k is the opacity, which, for simplicity, 
we take to be a constant, 0.03 cm 2 g" 1 . For our initial state 
we consider the atmosphere of pure helium with column depth 
y = 5.4 x 10 8 g cm" 2 at the point of ignition. The temperature 
at the start of the runaway (determined as the point where tem- 
perature derivatives of heating and cooling match) is 7g = 1 .64. 
The height of the top fluid layer, which in our model represents 
temperature, is normalized in units of the scaleheight of the un- 
perturbed atmosphere. In order to set off the runaway we raise 
the temperature at one location on our grid on a scale much 
smaller than the Rossby radius. The sequence of snapshots 
of subsequent evolution of layer height (temperature), instan- 
taneous net heating Q(x,t), and ageostrophic and geostrophic 
velocities in the top layer is shown in Fig. |] in the one-layer 
limit of our model (e = 0). 5 In this simulation, the Coriolis pa- 
rameter is constant everywhere on the grid and is representative 
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is thus set by the magnitude of the Coriolis force which turns 
the cross-front fluid motion into the flow parallel to the front. 
Since the scaleheight inside the front can change by factors of 
order 10 for strong bursts, nonlinear effects become important. 
Typically, however, the characteristic front thickness is of the 
order of the Rossby radius in the hot material behind the front 
[~ 3km for our parameters, cf. Eq. (^4j) ]. Since the burning 
must be complete inside the front width, the speed of the front 
should be such that it moves one Rossby radius in a character- 
istic nuclear time. For our parameters, it takes 0.1-0.2 seconds 
to go through the fast part of the 3a burning, which yields an 
estimate of 15-30 km s -1 for the front speed. Behind the loca- 
tion of peak heating the temperature continues to rise, but the 
magnitude of the pressure gradient starts to decrease, and the 
isallobaric wind is directed opposite to the motion of the front. 
The relative amount of cross-front fluid motion in the forward 
and backward directions can dramatically influence the speed of 
propagation of the front and depends on frictional effects due to 
(turbulent) viscosity and on drag due to momentum-conserving 
transport of fluid between layers. We shall refer to the latter 
effect as momentum coupling. 

Figure || shows the internal structure of two fronts com- 
puted in the one-layer approximation with identical heating 
functions but differing in the strength of momentum coupling. 
Front A (fig. ^a) is computed without the terms proportional 
to (v'2 AjV - V\x,y) m Eqs. (fiO])-(|43|) (momentum-nonconserving 
front), while front B (figTgb) is computed with momentum- 
coupling terms retained. Viscous friction is turned off in both 
cases (/x/ = 0). Physically, case A represents the situation where 
the fluid from the cold layer is being injected into the hot layer 
at the velocity of the hot layer, while in B the fluid is injected 
maintaining the velocity of the cold layer (which is zero for 
one-layer case e <C 1)- In the latter case, as the injected fluid is 
accelerated to the velocity of the hot layer, it exerts drag on the 
hot layer, which is reflected in the friction-like velocity depen- 
dence in the coupling terms in Eqs. (|40"|)-( 45 ). 6 

The two fronts move with different speeds: front A achieves 
20 km s" 1 , while front B has a speed of 60 km s" 1 . The 
main qualitative difference in the structure of these fronts is 
the variation of ageostrophic velocity v x with position inside 
the front. For front A the cross-front displacement of a fluid 
element is symmetric, with a fluid particle returning to the 
same x coordinate after the front passes, while for front B the 
ageostrophic component is asymmetric, skewed towards higher 
forward velocity. With stronger cross-front circulation front 
B achieves a larger speed of propagation. The cause of the 
larger ageostrophic velocity in case B is the modification of the 
geostrophic balance brought by the momentum-coupling drag. 
The isallobaric wind relationship is now v x = ~(§// 2 )[^V/! + 
{Q/h)Vh\. Since Q/h ~ l/t„ the two terms in this formula are 
of the same order of magnitude, and the cross-front wind is en- 
hanced where the magnitude of the pressure gradient is growing 
with time (ahead of the peak heating), and diminished where 
the pressure gradient is falling. As the heat transport depends 
on the cross-front circulation, the momentum-conserving front 
should exhibit larger speed of propagation. This is consistent 
with the results obtained in Sec. || where it was shown that, in 
general, friction can increase the front speed. 

The nature of the temperature overshoot in Fig. ||a can also 
be attributed to the differences in ageostrophic velocities be- 



tween the two fronts. In front A, the hot fluid flows towards the 
back of the front as much as it does in the forward direction. 
This flow increases the temperature of the fluid at the back of 
the front beyond the equilibrium value. After the helium fuel 
is depleted the material cools by radiation, and a "cooling tail" 
develops because of the time delay in the start of cooling due 
to the finite speed of the front. Momentum conserving fronts 
also develop cooling tails, but because of the larger propaga- 
tion speed the front moves through a larger distance before the 
fuel depletes e noug h for cooling to dominate heating (see dis- 
cussion in Sec. 4.3). 

We can derive a formal expression for the front propaga- 
tion speed by transforming equations (|40|-p3|) for the one-layer 
model to the frame comoving with the front. Expressing the 
ageostrophic component v x as a sum of constant front speed V/ 
and a residual v' x , and assuming that in the comoving frame the 
front is not evolving in time we obtain: 



G, 



d_ 

'dx 
d_ 

! dx 



"^- g ¥x h+fVy 

Vy=-f(v' x +Vf)- 



h 



(v' x +Vf), 



(65) 



(66) 



(67) 



As can be seen from Fig. pi the geostrophic velocity v y reaches 
an extremum inside the front. Denote all values at this point 
with an asterisk. At the extremal point, v' x \* = -vj -v*(Q{h*) + 
Hf)/h*, and we can solve the above system of equations for the 
speed of the front: 



Vf = - 



Q{h*)-h* 



Q(h*)+f-i f 



dx' 



f 



(68) 



Momentum coupling and viscous friction effects are contained 
in the second term on the right hand side of (|68|). Without 
them the speed of front propagation is controlled by the ef- 
fective heating and the pressure gradient at the peak of the 
geostrophic velocity. The effective heating Q(h*) - h*J^v x \* 
is the nuclear energy generation minus the heat carried by the 
cross-front flow. When the pressure gradient at the peak is ap- 
proximated as h* /or we get that the front speed is set by the 
characteristic front width divided by the effective burning time, 
as in (^|). The momentum-coupling drag term is of the same 
order of magnitude as the first term in (g8|), and increases the 
speed of propagation and the width of the front by a factor of 
2-3. 

3.6. Multilayer dynamics 

Although illustrative, the one-layer model represents the mo- 
tion of only the hot layer of fluid and does not capture all of 
the dynamics of the front. In order to better understand the 
propagation of the deflagration front through the neutron star 
atmosphere we numerically solve the time-dependent system 
of equations (|40|-p3|) for two layers of isentropic fluid. Fig- 
ure 1 1 demonstrates the structure of the developed momentum- 
conserving burning front without interlayer friction propagat- 
ing from left to right. In this particular run the density contrast 
e was set to 0.2, with the initial thickness of the lower layer 



6 Such drag converts some energy from the flow in the hot layer into heat. According to our model, heating represents injection of fluid into the top layer, so, if 
one is to be precise, the heating function should be renormalized to account for the frictional heating. When included in simulations, however, such renormalization 
accounted for no more than a 25% increase in front speeds, while not affecting the qualitative picture. Henceforth we will ignore the frictional heating. 
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chosen to be three times the thickness of the top layer. The 
first panel shows the evolution of layer thickness through the 
front. As the fluid is pumped to a higher entropy state, layer 2 
expands and layer 1 contracts; the ratio of contraction to expan- 
sion is proportional to the density contrast. We take the hori- 
zontal location of the peak in the heating rate as the center of 
the front. From the plot of the ageostrophic velocities we see 
that there is a divergence of the fluid at the top level and a con- 
vergence at the bottom level. The low-entropy fluid is drawn 
towards the center of the front and after ignition outflows in the 
high-entropy layer. The motion of fluid at the lower layer can 
also be interpreted as a response to the gradient in the column 
depth across the front created by divergence of the light fluid on 
top. In addition to creating convergence towards the center, the 
overpressure ahead of the front center also forces some flow to 
go in the direction of front motion. This "snowplow" effect will 
become important in the regime of strong top-bottom coupling; 
see the next subsection. 

The dominant component of motion in both layers is the 
geostrophic velocity (third panel of Fig. |]a). The direction of 
this flow is different between the layers: divergence at the top 
generates anticyclonic motion, while convergence at the bottom 
generates cyclonic motion. The fluid velocity in the upper level 
is larger than that in the lower level. The relative velocities in 
the two layers depend on the density contrast and on the inter- 
layer friction. 

3.7. Frictional effects 

We have argued in section || that friction between the top and 
the bottom of the atmosphere modifies the geostrophic balance 
within the front, and changes the velocity of the flow across the 
front. The cross-front velocity attains properties of a diffusive 
flow governed by the diffusion constant (plj). Of course, the 
diffusion analogy is only formal and refers to the spreading of 
pressure gradients with time in the presence of friction. De- 
pending on the strength of frictional coupling, the speed of the 
front may be either enhanced or diminished [cf. Eq. (|3~5|)]. Re- 
call that in our analytical estimates (§ we parametrized the 
strength of friction by the coupling time ff r , while in our two- 
layer model, described by eqs. (BQ-Ba), the friction strength is 
parameterized by fif/h, where [if is a constant. To facilitate 
the comparison between analytical and numerical results, we 
define, for our one-layer model, the characteristic timescale of 
frictional coupling as f& = h*/fj,f, where h* is the layer height 
at the peak of the geostrophic velocity inside the front (in the 
presence of friction this is approximately but not exactly the lo- 
cation of peak heating). In Fig. || we plot the speed of the front 
obtained from one-layer simulations versus the dimensionless 
friction parameter fi = (?f r /) _1 and show the fit using the ana- 
lytic formula for the front speed in Eq. (|34]). 

There are several distinct regimes depending on how the 
timescale of frictional forcing compares to the characteristic 
timescales of the problem. When l/f<t„< t& (or p < 0.01 
for our parameters) the effects of friction are insignificant and 
the speed is close to the value obtained in section 3.5 after in- 
cluding momentum coupling (~ 60 km s -1 ). For 1// < ff r < t„ 
(0.01 < p < 1) the speed increases as \/p/(l + p 2 ) in agree- 
ment with (f34]). The friction in this regime acts to reduce the 
geostrophic velocity, which extends the front over a larger spa- 
tial scale. The cross-front diffusion constant ( |oT| ) is maximum 
when friction acts on the rotation timescale p = 1, and at this 
value the front speed reaches the peak value (464 km s~ for 



our parameters). The structure of the two-layer model at p = 1 
is shown in Fig. ||b. The main distinction of this case from the 
case of zero friction (fig. ||a) is the form and the larger magni- 
tude of the ageostrophic speed, which implies very strong cross- 
front circulation. In the top layer the flow is purely in the for- 
ward direction with a strong return current in the bottom layer. 
The value of the geostrophic speed is clearly reduced compared 
to the no-friction case, and the strong coupling between the lay- 
ers causes the flow on the bottom to have the same anticyclonic 
direction as the flow in the top layer. 

When the friction is further increased (pi > 1), the cross-front 
diffusion ( |6r} ) diminishes and the front speed decreases. The 
cross-layer coupling now acts to reduce the ageostrophic com- 
ponent of velocity by making both layers move together. This 
"snowplow" effect eventually stalls the front when the friction 
is very large. Rather than create circulation across the front, 
the pressure gradient between the hot and cold material tries 
to push forward the whole cold atmosphere, which makes the 
front steepen up and slow down. 

It is interesting to speculate whether the extraordinarily high 
front speeds of hundreds of km s" 1 observed in simulations can 
be expected to occur in reality. One does not expect a well- 
defined burning front to arise if f prop < t n , where f prop ~ 7r^/vf ront 
is the timescale on which the formed front would cross the half- 
circumference of the neutron star. In other words, it makes no 
sense to talk about a well-defined front with a width A ~ Vf ront ?„ 
that is larger than the size of the star. This implies that if 
Vfront ^ iTR/t n ~ 300 km s _I , the front never forms and the star 
is essentially ignited simultaneously. Depending on the physi- 
cal source of friction in the neutron star atmosphere there could 
then be two possible scenarios for the spread of nuclear burning. 
If the strong friction is temperature dependent, and increases 
when the local perturbation grows into runaway (as may be the 
case for convective friction), the runaway will begin with no 
friction present and, therefore, will be confined to a hot spot 
of the size of the Rossby radius. As the runaway progresses, 
the friction will increase and the temperature perturbation will 
quickly (~ 0.1 s) spread over the neutron star surface. If, how- 
ever, the friction is strong and present all of the time during 
the runaway (e.g., if a magnetic field is threading the fuel), an 
initial hot spot may not appear at all. Rather, the temperature 
perturbation spreads away from the heat source faster than the 
local thermonuclear runaway can develop, and the whole sur- 
face is likely to ignite almost simultaneously in a spherically 
symmetric fashion. 

4. GLOBAL HYDRODYNAMICAL FLOWS DURING X-RAY BURSTS 
AND CONNECTION TO OBSERVATIONS 

4.1. Likely location of burst ignition 

So far we have assumed that the pre-burst conditions are 
identical everywhere on the neutron star. However, because the 
star rotates, the effective gravity felt by the fluid elements near 
the equator is somewhat (up to ~ 25%) smaller than that felt 
by the fluid elements near the poles. As we now argue, this 
asymmetry implies that, even if accretion is perfectly spher- 
ically symmetric, the fuel near the equator is likely to reach 
ignition conditions first. 

We assume that, prior to a burst, the accreted material is 
brought into corotation with the rest of the star either by effec- 
tive hydrodynamic viscosity produced by Rayleigh-Taylor or 
baroclinic instability (Fujimoto 1988), by weak hydromagnetic 
stresses, or even by microscopic viscosity (Cumming and Bild- 
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sten 2000). The timescales for the first two processes are small 
compared to the interval between the bursts, so this assumption 
is probably accurate as a zeroth-order approximation. We also 
assume that the gas deposited onto the star can redistribute to 
achieve hydrostatic balance 7 . 

In hydrostatic equilibrium, the pressure at the bottom of the 
accreted ocean should be the same everywhere, 



Pbottom = m(A)geff(A) = const, 



(69) 



where m(X) is the accreted column density, \ = ir/2-8 is the 
latitude, and g e ff = |V^| is the effective gravitational accelera- 
tion. Here cf> is the sum of the gravitational potential, cf) gI , and 
the centrifugal potential, SI 2 /? 2 cos 2 A/2. The effective gravity 
geft on the equator is less than that on the poles due to rotation 
of the star: <5g e ff/geff ~ (^/^k) 2 , where u>k is the Keplerian an- 
gular frequency at the neutron-star surface. For a neutron star 
rotating at 300-600 Hz, the relative difference is a few percent 
(up to 1/4). Therefore, the column depth m(X) to a surface of 
a given pressure is a factor of <5g e ff/geff higher on the equator 
than on the poles. Now concentrate on two fluid elements, both 
at a pressure /bottom, one of which is located on the equator 
while the other one is somewhere at a latitude A ^ 0. Imagine 
that an additional amount of gas is accreted, and is allowed to 
come to corotation and redistribute itself on the surface of the 
star in order to achieve hydrostatic balance. Our fluid elements 
will now be compressed to a pressure p\ >otlom + Ap (since they 
are assumed to be in hydrostatic balance, they must have the 
same pressure). But, as evident from Eq. (|9|), the increase in 
the column depth, Am(A), must be larger on the equator than 
away from it. Therefore, fuel arriving on the equator will reach 
a given column depth faster than the fuel located off-equator. In 
other words, the local accretion rate is inversely proportional to 
the local effective gravity strength, 



m(X) = rho 



go 



geff(A) ' 



(70) 



where = M/ (AnR 2 ) is the average local accretion rate, and go 
is the average acceleration of gravity, chosen so that the integral 



over the surface of the star J^ 2 ^ 1/2) cos X[go/g(X)]dX 
The local effective accretion rate seen on the equator is at least 
a few percent higher than that seen at the poles. 

At a given local accretion rate, which determines the thermal 
profile of the accreted ocean, a certan critical column depth is 
necessary for the thermal instability and a thermonuclear run- 
away (see Bildsten 1998 for a review of spherically symmetric 
ignition conditions). As we argued above, the critical column 
depth necessary for a runaway will be achieved first close to the 
equator. We can therefore expect that ignition happens close to 
the rotational equator. 

4.2. Burning front propagation on a (3-plane 
The speed of the front propagation is larger near the equa- 



tor, where / — > 0, than near the poles, where / — * f m 



(cf. Eqs. [ |2"4] ] and [34]]). To understand qualitatively the ef- 
fect of the latitude dependence of the Coriolis parameter, we 
model front propagation on an x — y plane oriented with x axis 
east, and y axis - north. The Coriolis parameter is taken to be 
/ = (3y, where = 2il/R. This approximates a rotating sphere 
in the equatorial region, and is known in geophysical literature 
as the equatorial /3-plane. 

1-D simulations of equatorial ignition. We begin in 1-D 
with the simulation axis oriented along the y axis from the 
equator to the north pole (and allow variation in quantities only 
along this direction). We solve a one-layer shallow water model 
without friction, which is initialized by increasing the temper- 
ature on the equator. A series of snapshots in time is shown in 
Fig. ^. Since the Coriolis parameter is zero on the equator, it 
is not possible to establish geostrophic balance there. Gravity 
waves from the initial perturbation then propagate towards the 
pole until they are reflected from the region with a finite Cori- 
olis parameter at approximately y = ore = (Vgh/P) 1 ^ 2 - This 
distance is called the equatorial Rossby radius and is the char- 
acteristic width of the equatorial waveguide for gravity waves. 
The trapped gravity waves are amplified by thermonuclear en- 
ergy release on each pass. For our choice of initial conditions, 
after about t n "JgH/a.RE = tniP^/gh) 1 ^ 2 ~ 10 2-3 passes the whole 
equatorial waveguide region undergoes a runaway and two 
geostrophically supported burning fronts propagate towards the 
poles as "walls of fire". The propagating fronts steepen and 
slow down as expected because of the increasing /. For burn- 
ing on a rotating neutron star the width and speed of the front 
decrease by a factor of 2^l/{(3y/gh) l l 2 w 4 as the front propa- 
gates from the equator to the pole. 



2-D model for front propagation. 

infer the front speed on the /3-plane: 

VHam e (x,y) = V e 



From Eq. (B4j) we can 



(71) 



v2 



where x = ((3tf r )x, y = ftf r = ((3t&)y, and v eq is the speed of 
the front at the equator. We have assumed for simplicity that 
tft < f„; if this is not the case, then in the above expressions f n - 
must be replaced by ff r f„/(f„ + r/tfX 



As was discussed in Section 3.3, ignition happens over a 



patch of a finite size which is determined by the Coriolis param- 
eter and the strength of friction in the atmosphere. However, for 
the current discussion, we assume that ignition is highly local- 
ized near some point on the /3-plane. In Fig. || we show how 
the front line develops depending on the location of the igni- 
tion point, yo = /off r , where fo is the Coriolis parameter at the 
ignition location. To compute the evolution of the front line we 
used the Fast Marching Method on a triangulated mesh (Sethian 
and Vladimirsky, 2000) with speed function given by (J7l|). 9 

Below, we comment separately on the case where the ignition 
point is on the equator, and on the case when the atmosphere is 
ignited at some non-zero latitude. 



7 Inogamov and Sunyaev (1999) have considered the spreading of material accreted onto a rotating neutron star from a thin equatorial disk. They have concentrated 
on the hot radiation-pressure dominated fluid at the very top of the atmosphere, which has not had time to frictionally couple to the rest of the star. The nuclear fuel 
participating in an X-ray burst is below the hot radiating layer considered by Inogamov and Sunyaev. We therefore expect that the mid-lattitude "rings of fire" found 
by Inogamov and Sunyaev will not affect the location of burst ignition. 

8 Our equations break down at the equator, because the horizontal Coriolis force due to vertical motion, 2ncosAi>; fs 2f!cos \{h/t n ), is no longer negligible 
compared to the Coriolis force due to horizontal motion, fv g f» 2f! sin A(g/i)'/ 2 . However, this breakdown occurs only for latitudes A < (h/g) [ / 2 (l/t„) ~ 1(T 5 (for 
typical values of h and f„), i.e., very near the equator. 

9 We thank Alexander Vladimirsky for introducing us to this powerful method and helping to produce Fig. 
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Equatorial ignition [Fig. ||aJ7- The front is roughly a cir- 
cle until it reaches \y\ ~ 1, i.e., until it reaches the latitude 
where / w 1/ff, and the front speed is significantly reduced. 
Burning then quickly spreads along the equator, and the front 
line is nearly parallel to the equator when the flame reaches 
9 = ftfi ~ 6. For a typical strongly accreting neutron star in an 
LMXB the Coriolis parameter at the pole is 4-7 x 10 3 rad s" 1 , 
and the typical frictional coupling strength l/ff r is anywhere 
between 10 s -1 for pure momentum coupling, and 10 3 s" 1 for 
frictional coupling by strong convective turbulence. Our simple 
2-D model then indicates that the front becomes parallel to the 
equator as soon as it reaches the mid-latitudes, or perhaps even 
closer to the equator. 10 

What is the observational significance of this result? The 
asymmetry between north-south and east-west front propaga- 
tion speeds, and the dependence of the speed on the convective 
coupling strength may explain why only a fraction of X-ray 
bursts shows detectable nearly coherent oscillations during the 
burst rise. Indeed, during the burst rise, while the flame, ignited 
at a spot, is still spreading around the star, rotational modulation 
of the X-ray flux should be evident. Yet it is observed only in 
some bursts. A possible explanation previously proposed in the 
literature (Miller 2000 and references therein) is that the oscil- 
lations in the X-ray flux are washed out when the burst X-rays 
scatter off a hot cloud which is thought to surround strongly 
accreting neutron stars. This explanation, however, does not 
resolve why there are oscillations in some bursts and no oscil- 
lations in others. Muno et. al. (2001) recently concluded that, 
in fast (~ 600Hz) rotators (as inferred from the burst oscillation 
frequency 11 ), only the strong radius-expansion bursts show os- 
cillations during the rise, whereas for slower rotators (~ 300Hz) 
there is no such correlation. 

In our model, the equatorial front speed is larger than that at 
the poles, and, as argued in § X. 1 , the burst is likely to ignite 
at the equator. If the disparity between the front speed at the 
equator and the poles is large (i.e., if 2f2ff r ^> 1) then the entire 
equatorial belt is likely to ignite in the beginning of the burst, 
and the asymmetry in burning which is needed for burst oscil- 
lations may disappear. Increasing the burst strength enhances 
the frictional coupling (makes ff r shorter) and makes the burn- 
ing pattern more asymmetric, thus making the burst oscillations 
during the rise more likely. This means that equatorially ig- 
nited strong bursts are more likely to be asymmetric than weak 
bursts. What about the dichotomy between fast (600 Hz) and 
"slow" (300 Hz) rotators? For slow rotators, the red uction in 
equatorial gravity due to rotation, discussed in § 4. 1 , is a fac- 



tor of 4 smaller, so there may be a significant probability that 
bursts ignite off the equator, and, hence, are more asymmetric 
regardless of burst strength. 

Nonequatorial ignition [Fig)Q (b), (c)]. The front is roughly 

G2 



4.3. Zonal flows and frequency drifts 

A salient feature of burst oscillations is their presence in 
the tails of bursts, presumably after the entire ocean has been 
burned, as well as the increase of the oscillation frequency as 
the star cools. We speculate on the possible origin of the inho 



mogeneity in the ashes that gives rise to the oscillations in § \A 



here we discuss the implications of our global X-ray burst sce- 
nario on the frequency of these oscillations. In addition, we 
predict that yet unidentified oscillations may be present during 
the burst rise as well. 

The burning front leaves hot ashes in its wake. These hot 
ashes then cool on a characteristic timescale f 000 i, defined as 
dhhot/dt ~ /ihot/fcooi- If burning starts in the equatorial region, 
as we have argued, then the equator has a headstart in cooling, 
and in the cooling wake the equatorial temperature will be the 
lowest, increasing towards the poles. This temperature gradi- 
ent will drive a zonal thermal wind directed backwards relative 
to the neutron-star rotation. If an inhomogeneous feature (e.g., 
a vortex-see the next subsection) is trapped in this backward 
zonal flow, the frequency of the flux modulation due to this fea- 
ture will appear lower than the neutron-star spin frequency. The 
speed of the backward flow can be estimated as follows: 
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The drift frequency corresponding to the flow is given by 
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In the mid- latitude Vfl ame ~ /?/f r ise, where f r i se is the timescale of 
the burst rise. The angular speed of the drift in the mid-latitude 
is then 
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here Vk and v s are surface Keplerian and spin frequency, respec- 
tively. For comparison, the frequency due to the radial expan- 
sion of the burning layer (Strohmayer et. al. 1998, Cumming 
and Bildsten 2000) is 
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The geostrophic drift from Eq. ( [75] ) and the lift-up drift from 
Eq. ( f75[ ) are independent of each other, since the lift-up drift is 
related to the horizontal component of the spin vector which we 
have neglected in our calculations. Both drifts must be present 
in real bursts. One effect is not obviously dominant over the 
other; the "free" parameters in Eq. ([75]) decide which effect is 
more important for a particular burst. 

Recently Van Straaten et. al. (2001), Galloway et. al. (2001), 
Wijnands et. al. (2001), and Cumming et. al. (2001) showed that 
the lift-up drift is insufficient to explain the observed magnitude 
of the frequency drifts in some bursts. Perhaps the geostrophic 

10 When the front is parallel to the equator, the flows associated with the front are also parallel to the equator; such flows are referred to as zonal currents in 
geophysical literature. In this_naper we do not consider in detail the nonaxisymmetric instabilities which may be present when the strength of the zonal current is 
latitude-dependent, but see § K4 

"For one of the bursters, there is evidence that it spins at half of the burst oscillation frequency (Miller, 1999; see also Strohmayer 2000). 



a circle so long as the distance to the ignition point \J Ax 2 + Ay 
is less than 1 . Beyond this distance, the front becomes increas- 
ingly deformed. As the edge of the flame approaches the equa- 
tor, the front accelerates, and the equatorial belt is quickly ig- 
nited. After this, the evolution is similar to the equatorial igni- 
tion case discussed above. 
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wind could provide the missing piece. A significant amount of 
modeling of burst light curves is required to make a more de- 
tailed comparison with observations; this is the subject of our 
current research. 

Note that from Eq. ([73]) one can infer significantly larger 
drifts close to the pole and the equator, than the mid-latitude 
drifts given by Eq. (75). However, close to the equator the 
geostrophic balance is broken, and Eq. ([73|) is no longer valid. 
The origin of divergence of Wdrift at the pole is that we formally 
take the distance to the pole to zero, while keeping non-zero 
Vflow Perhaps, physically one should not consider distances to 
the pole closer than the width of the burning front. 

Like the lift-up drift, the geostrophic drift velocity asymp- 
totes to zero as cooling continues, since the scaleheights at the 
equator and the pole equalize. Thus, the final X-ray modula- 
tion frequency asymptotes to the neutron-star spin. A possible 
observational signature that can distinguish the two effects is 
the time dependence of the drift frequency. Consider an atmo- 
sphere undergoing black-body cooling in the tail of the burst. 
The scaleheight of the atmosphere then decreases with time as 
h(t) = /ihot(l +t /t c )~ 1 ^, where t c = c p ny 2 j '(3acr o 3 ) is the char- 
acteristic cooling time at the maximum expansion when the 
temperature is 7J). For Tq = 10 9 K, y = 5.4 x 10 8 g cm" 2 , and 
k = 0.03 cm 2 g" 1 , this cooling time is 2 s. Since the geostrophic 
speed depends on the gradient of the scaleheight, the X-ray os- 
cillation frequency due to geostrophic drift should change with 
time as: 
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where w^-m i s tne drift frequency at the start of cooling. Here 
we assumed that the feature that gives rise to the X-ray oscilla- 
tions remains at the same latitude for the entire duration of the 
cooling tail of t he b urst and is carried by the geostrophic flow 
(however, see § 4.4 regarding the possibility of latitudinal drifts 
which can modify Eq. [f76|]). 

On the other hand, the frequency drift due to the radial expan- 
sion of the atmosphere depends on the scaleheight of the col- 
umn, rather than its gradient. The oscillation frequency asymp- 
totes to the neutron star rotation frequency as 1— (u;{j ft _ up /fi)[l + 

t/tcY 1 / 3 - Since it is likely that in real bursts both the lift-up and 
the geostrophic drifts are present, the X-ray oscillation signal 
will have a time-dependence corresponding to a superposition 
of the two effects. If both drifts are of the same order of mag- 
nitude, different time-dependence of the two effects will result 
in the initial stage of the frequency evolution being dominated 
by the geostrophic drift, while the late stage will be dominated 
by the lift-up drift. We note that the time scalings given above 
assume a very simplistic model of the evolution of the atmo- 
spheric scaleheight during the cooling. Clearly, more detailed 
simulations are necessary for robust comparisons with observa- 
tions. 

As demonstrated in Sections |]and|3|, the dominant flow dur- 
ing the burst rise is a geostrophic current directed parallel to 
the front line. For the moment, let us assume that the friction 
within the front is weak. Then the geostrophic current is given 
by v g ~ (ghhot) 1 ' 2 ~ 1.5-3 x 10 8 cm s" 1 , and is independent 
of latitude [this follows from eqs. ( |l9[ ) and If there is 

a temperature inhomogeneity entrained in this current (similar 
to the inhomogeneities responsible for oscillations in cooling 
tails of bursts which were discussed above), we may expect a 
flux modulation at a frequency higher than the spin frequency; 
it will appear as an upper sideband of the main burst oscillation 



frequency during the rise. When the front is at a latitude A, the 
frequency of this sideband is v^/(27rl?cosA). As the flame prop- 
agates towards higher latitudes and carries the inhomogeneity 
along, the frequency chirps up. 

The time evolution of this chirp is easy to estimate. The 
speed of the front depends on latitude as v/ = v^ ' 6 / sin A, where 

Vy° le ~ y/ghhot/ (2£lt n ) is the speed of the front near the pole. 
Expressing the front speed as v/ = RdX/dt, we solve this equa- 
tion for the front latitude as a function of time: cos A = l—t/t p , 
where t p = R/v p ^ e . The time evolution of the sideband fre- 
quency is then 
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Here, 



chiip 



(ghy/ 2 /2nR is the sideband frequency in the re- 
gion of equatorial belt, which, for typical parameters in our sim- 
ulations, ranges from 25 to 50 Hz. The final frequency of the 
chirp depends on how close to the pole the front stops propa- 
gating. If we take this halting distance to be one Rossby radius 
fl P ole r*j 1 km, the final frequency of the chirp is ten times the 
equatorial value, or 250-500 Hz. 

When friction is present, the speed of the zonal flow within 
the front is 

^geostrophic ^ r a ^ r J \ ' 

J A /fnVflame 



where Vfl ame is given by Eq. (34). We have then 
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Close to the equator, where the Coriolis parameter / is small 
compared to the rate of frictional coupling 1 / ff r , the chirp fre- 
quency decreases as the front moves away from the equator. 
However, as soon as the front reaches the region with / ~ 1 / ff r , 
the chirp frequency increases as the front moves further to- 
wards the pole. When / is comparable to or greater than l/ff r , 
the magnitude of the chirp frequency is smaller by a factor of 
~ V ' tn/tfr than what it would be if the friction was absent. If 
measured, the time-dependence and the magnitude of the chirp 
frequency could be used to discern the importance of the effec- 
tive viscosity during the burst. 

4.4. Vortices and oscillations in the tail 

So far we have not addressed the nature of the perturbations 
that lead to the oscillations in the X-ray flux in the burst tail. 
The speed of the cooling flow described in the previous sub- 
section is latitude-dependent [cf. Eq. (j72|)]. Thus we can ex- 
pect a strong zonal shear, of the type that is observed in the 
atmospheres of giant planets. Such a zonal shear is known to 
be unstable to formation of vortices; sometimes — like in the 
case of Jupiter — these vortices become large and occupy a 
significant fraction of the atmospheric surface. We speculate 
that these vortices do form in a cooling wake of an X-ray burst, 
and are responsible for modulation of the X-ray flux in a burst 
tail. Currently we are constructing a 2-D shallow water simula- 
tion to address this scenario. Two-dimensional simulations will 
also allow us to study several other important effects which are 
not contained in the 1-D model. In particular, in the presence 
of a meridionally varying Coriolis parameter, vortices tend to 
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acquire a drift even if there is no background flow. This ef- 
fect is known as "/3-drift" in the geophysical literature (see, e.g. 
Chan and Williams, 1987), and is responsible for the north- 
west drift of tropical hurricanes on the Earth. The direction of 
the drift depends on the sense of rotation of the vortex, with 
cyclones drifting northwest and anticyclones southwest. In or- 
der to estimate the velocity of this drift we use the empirical 
formula from Smith, Li and Wang (1997): vp-^ ~ r m \/v m (5, 
where r m is the radius where the maximum fluid velocity v m 
relative to the vortex guiding center is reached. For a "hot spot" 
(anticyclonic) vortex with a radius equal to the Rossby radius 
(~ 1 km), and maximum internal speed v m ~ 10 8 cm s -1 on 
a star with (3 = 4 x 10~ 3 cm" 1 s" 1 , the drift speed is of order 
600 km s -1 . Depending on the direction of drift on the surface 
of the star this effect can yield an X-ray oscilation frequency 
of up to 10 Hz lower than the spin frequency of the star. It 
remains to be seen how burning inside the vortex and potential 
nonaxisymmetric instabilities (e.g. development of spiral arms) 
modify this estimate and affect the evolution of vortices in two 
dimensions. 

4.5. MHD dynamo and coherence of pulsations 

The burning front may be an ideal environment for an MHD 
dynamo 12 : it is likely to have both turbulence and, at least ini- 
tially, strong shear; additionally, the typical edde overturn time 
is much less than t„. The equipartition value of the magnetic 
field is 
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(80) 

where, for concreteness, we have assumed that the turbulence 
is of convective origin. The coherence length of the dynamo- 
generated magntetic field is of the same order of magnitude as 
the ocean scaleheight. 

As discussed in the Introduction, oscillations observed in 
the tails of X-ray bursts are highly coherent, with Q's of few 
thousand over the duration of the burst (van der Klis 2000). 
Both lift-up (Strohmayer 1998; Cumming and Bildsten 2000) 
and geostrophic models for the drift have difficulty accounting 
for this coherence: the burning ocean is strongly sheared, with 
the top of the scaleheight moving at a different speed than the 
bottom. Cumming and Bildsten (2000) argued that convection 
might enforce the vertically rigid rotation of the burning ocean, 
although they pointed out that this is far from certain since con- 
vective turnover time is comparable to the spin period. Numer- 
ical models of bursts indicate that strong convection is indeed 
present when the ocean is ignited and rises, but dies quickly 
when the fuel is exhausted and cooling begins. It therefore 
seems unlikely that the coherence of burst oscillations in the 
cooling tail of the burst could be accounted for by convectively 
enforced vertically rigid rotation. 

We propose a different idea. The B-field generated by the 
burning front will dynamically couple the top and the bottom 
of the cooling ocean 13 on the timescale (see, e.g., section 4 of 
Cumming and Bildsten 2000) 
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12 We thank Maxim Lyutikov for alerting us to this fact. 

13 It seems reasonable to assume that the small-scale magnetic field, with coherence length of 
the large-scale (~1 km, roughly the width of the front) lateral shear. 



which is typically much shorter than the duration of the os- 
cillations. Moreover, the coherence of the oscillations, Q ~ 
fcooiAcoupie is of the order of magnitude of observed values. We 
therefore argue that the vertically rigid rotation necessary for 
the burst coherence may be enforced by the small-scale mag- 
netic field generated by the dynamo in the burning front sweep- 
ing through the ocean during the burst rise. Note that after 
the burst this small-scale magnetic field is confined to the new 
ashes, while the freshly accreted matter can remain unmagne- 
tized. 

5. CONCLUSIONS 

We have constructed a new model for the spreading of defla- 
gration fronts during type I X-ray bursts. In contrast to previous 
models, we take into account the horizontal hydrodynamical 
flows arising due to the radial expansion of the burning atmo- 
sphere/ocean and the action of the Coriolis force due to rapid 
rotation typical for strongly accreting neutron stars in LMXBs. 
Our mechanism of heat transport relies on coherent hydrody- 
namical flows across the front that are set up as the fluid at- 
tains momentum balance, and achieves speeds of front propa- 
gation in excess of tens of kilometers per second, necessary to 
account for the observed sub-second rise times of X-ray bursts. 
Previous workers invoked large-scale convective turbulence (on 
scales much larger than the vertical scaleheight) in order to ob- 
tain comparable front speeds. The speed of the front is analyt- 
ically estimated in Eq. (^4|) for the case when convective vis- 
cosity is not important, and in Eq. (Q) for the case when the 
layers of the burning ocean are frictionally coupled (by, e.g., 
convection). 

In addition to the analytical arguments, we constructed and 
numerically evolved a two-layer shallow-water model of the 
burning layers. Our simulations agree very well with our an- 
alytical estimates, and we found that the physics of propagation 
of geostrophically supported fronts is independent of the details 
of the heating function, as long as it is a single-peaked function 
of temperature. This agreement gives us confidence that our 
results are valid for models with more realistic microphysics. 

We have outlined the behavior of global hydrodynamical 
flows in the neutron-star atmosphere/ocean during the burst, 
and showed that these flows may explain many of the features 
of observed bursts. The very short rise times of X-ray bursts 
are easily accounted for if the speed of the burning front is set 
by the requirement of geostrophic balance. The lack of burst 
oscillations in many bursts may be due to the fact that the speed 
of the burning front is very dependent on the location on the 
neutron star surface, and favors rapid propagation along near- 
equatorial latitude bands that wipes out asymmetry, rather than 
simple spreading of a hot spot that maintains asymmetry. The 
rather large frequency drifts of burst oscillations in tails of some 
bursts can be accounted for only if the lift-up drift due to the ra- 
dial expansion (Strohmayer 1997, Cumming and Bildsten 2000, 
Cumming et al. 2001) is combined with the geostrophic drift 
due to zonal flows in the cooling wake of the burning front. 
We argue that the burning front will generate strong (~ 10 9 G) 
small-scale magnetic field, and that this magnetic field will en- 
force a vertically rigid flow of the ocean in the wake of the 



10 m (the ocean scaleheight) will not have a substantial effect on 
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burning front. The vertical rigidity of the flow explains the ob- 
served high degree of coherence of the burst oscillations during 
the cooling tails of bursts, after convection has subsided and 
can no longer account for the coherence. We conjecture that 
strong zonal currents during the burst may lead to the forma- 
tion of vortices of the type observed in the atmospheres of giant 
planets. These vortices may be responsible for X-ray flux oscil- 
lations in the burst tail. In addition, we predict the presence of 
yet unobserved chirps during rises of bursts. Detection of such 
chirps (which could have been missed in previous observations 
because of the large frequency spans that they cover) will tell 
us about the details of burning front propagation, such as its ve- 
locity or the latitudinal extent of the surface of the star that is 
covered by nuclear burning. 

Our model is incomplete in several respects. We showed 
that the burning front on a rotating neutron star is a site of 
strong vertical and horizontal differential shear, and that this 
shear flow can transport entropy across the front. However, we 
omitted the small-scale mixing and heat diffusion, and thus we 
have not showed that this heat can indeed be delivered into the 
cold fuel, heat it up, and ignite it as the burning front propa- 
gates. Since flows with strong shear tend to have both local and 
global hydrodynamical instabilities, we are confident that such 
heat exchange does occur and is robust. Only 2- or 3-D hydro- 
dynamical simulations will be able to authoritatively settle this 
issue, and we have outlined how our results would be modified 
if the mixing exists, but is not as efficient as we have assumed. 
In addition, such simulations will be useful for understanding 
the dynamical (frictional) coupling between the layers of the 



burning ocean, which, as we showed, has a substantial effect 
on the speed of propagation of the burning front. Finally, fur- 
ther modeling of the global nonaxisymmetric instabilities aris- 
ing in the horizontal shear flows may confirm our conjecture of 
trapped vortices as the origin of flux modulation at late times 
during X-ray bursts. 

In this paper we have demonstrated the usefulness of the 
shallow-water model for understanding X-ray bursts. Two- 
dimensional simulations utilizing the shallow-water model on 
a sphere will complement future vertically-resolved hydrody- 
namical simulations, and will allow us to generate realistic 
lightcurves and make more refined predictions of frequency be- 
havior of burst oscillations. 
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APPENDIX 



FORMAL DERIVATION OF MODEL EQUATIONS 

The shallow-water equations (|40|-p3|) for incompressible fluid can be derived as a formal limit of the Euler equations for compress- 
ible fluid. This is useful for establishing the connection between the entropy sources in compressible fluid and mass sources which 
are used to simulate heating in the incompressible case. Using the continuity equation and the ideal gas law, the heat equation i 
can be rewritten as: 

(1-k)— = -kT\7-v+—, (Al) 
at c p 

where k = 1 - 1 j^. In the limit of an incompressible fluid, 7 — > 00 and n — ► 1 . We can then equate the right hand side of (|A1[) to zero 



We consider a layer of uniform temperature Tq and height h, which in our model is equal to the scaleheight c p To/g. Integrating (Al) 
in the vertical direction over the layer height, we get: 



hV\ 1 • V+ h(v z \ z= h - v z \ z= o) ■ 



J 



(A2) 



We assume a hard surface at the bottom of the layer (v z |~o = 0), and use v z \ z= h = dh/dt as the velocity of the surface. This yields 
the continuity equation as in (^0[^lj) with effective mass source Q = J/g representing heating. The momentum equations ( p3[fi3| ) 
directly follow from the Euler equations for a constant density layer in hydrostatic balance. 
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FIG. 1 . — Illustration of a burning front moving from the hot to the cold region in the atmosphere/ocean 
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FIG. 2. — Front speed as a function of friction strengl 
the analytical expression for the front speed (equation ( 
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f-'h /y} ■ Dots represent the results of the shallow-water simulation (§p.7p and the solid line is the fit using 
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FIG. 3. — Adjustment to equilibrium in a fluid with an entropy gradient, a) Initial configuration with two columns of equal mass with uniform potential temperature 
6l > 02 in vertical hydrostatic balance. Upper surfaces are at the same external pressure p t , but not at the same height; b) Final configuration - entropy stratification. 
Higher entropy fluid is on top, and the interface is at pressure j(p, + p;,). 
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FIG. 4. — One-layer model of the evolution of a localized source of heating into two propagating burning fronts. Shown (from top to bottom): temperature 
(scaleheight) of the layer in units of 10 s K; instantaneous heating rate (in units of 5 X 10' 9 erg g -1 s -1 ); cross-front ageostrophic velocity v x and geostrophic velocity 
v y parallel to the front (velocities are in units of gravity wave speed of the cold material, 1.5 X 10 8 cm s -1 ). Time is increasing left to right with frames separated by 
0.125 sec. 
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FIG. 5. — Internal structure of burning fronts in the one-layer model. Shown are the temperature of the layer, the instantaneous heating rate (in units of 
5 X 10 19 ergg~'s~'), cross-front ageostrophic velocity v x , and tangential geostrophic velocity v y (velocities in units of gravity wave speed in the cold material, 
1.5 X 10 8 cm s -1 ). a) front without the drag due to momentum coupling; b) momentum-coupling drag included. 
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FIG. 6. — Internal structure of burning fronts in the two-layer model with e = 0.2. Shown are the thickness (temperature) of the layers with layer 2 (solid) added 
to layer 1 (dashed), the instantaneous heating rate (in units of 5 X 10" erg g _1 s - '), cross-front ageostrophic velocity v v , and tangential geostrophic velocity v y for 
two layers (velocities are in units of gravity wave speed of the cold material 1.5 X 10 8 cm s~'). a) Momentum-conserving front without viscous friction (fif = 0). 
For demonstration purposes, v x and v y for the bottom layer are increased by factors of 2 and 5 respectively; b) Momentum-conserving front with friction (fif = 1). 
v y for the bottom layer is increased by factor of 5. 




FIG. 7. — Ignition and propagation on /3-plane. Frames are separated by 0.06 s 
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